% This code can be executed to generate Figure 4A-C. This will plot the raw
% data (load), elastic (EF), and peeling (y1 and y2) force models for three 
% representative trials at 25 (C), 42.5 (A) and 57.5 (B) degrees. The 
% peeling maxima are purple dots (peelmax). The global maxima are dots 
% colored acording to their buckled conformation and junction angle (peakload).
% The point where the elastic and peeling energies are equal is marked with
% a yellow triangle (intersectUFelas). The point where the elastic model
% becomes linear is marked with an orange square (linearpts). The peeling
% maxima are also plotted as a function of the junction angle as a purple line
% (peelmax). Additional formatting is done externally to exactly match the
% published Figure 4.
% The folder titled, "Raw Data Test Files," must be open to run this code! 
clear;
%Initiate Tape Constants
t=0.058; %mm
w=12.7; %mm
Area=w.*t;
E=467.9745; % Pa
gamma=10.3014/1000;  %N/m
beta=mean([189.393 190.444 184.519 182.807 186.996]); %mm
TS_tape=23.2416;   %N
slope=(E*w)./1000; %N/mm
plotnorm=w;

%Determine which trials to analyze & plot
names=['2020-04-15_020_57.5.txt',...
    '2020-04-15_014_42.5.txt',...
    '2020-04-15_007_25.0.txt'];

%Offsets are needed to potentially remove slack from raw data to align with
%theoretical curves. This is done as guess & check.
dat_offset=[0.1 0.0 0.14 0];
text_offset=[0.5 0.5 -1];

%short algorithm to determine where each file name starts and ends to
%determine the junction angle information which is stored in the filenames
for i=1:length(names)
    chunk=i+3;
    test=names(i:chunk) == '.txt';
    if all(test)
        numchar=chunk;  % number of characters in each file name
        break
    end
end

%Initiate Variables
ntrials=length(names)/numchar;  % must be uniform to differentiate trials  40-41
jj=zeros(1,ntrials);
jj(1)=1;
kk=zeros(1,ntrials);
kk(1)=numchar;
angles=zeros(1,ntrials);
peakloads=zeros(10,ntrials);
peakload=zeros(1,ntrials);
xmaxs=zeros(1,ntrials);
pklocs=zeros(1,ntrials);
peelmax=zeros(1,ntrials);
linearpts=zeros(1,ntrials);
intersectUx=zeros(1,ntrials);
intersectUFelas=zeros(1,ntrials);
Lo=zeros(1,ntrials);
adjustloads=[];
adjustloc=[];

%Predetermine Each Trial's Color
colors=zeros(ntrials,2,3);
colors(1,1,:)=[70/255 0/255 130/255]; %data 1
colors(1,2,:)=[170/255 70/255 255/255]; %model 1
colors(2,1,:)=[0/255 120/255 80/255]; %data 2
colors(2,2,:)=[60/255 240/255 165/255]; %model 2
colors(3,1,:)=[70/255 150/255 180/255]; %data 3
colors(3,2,:)=[10/255 200/255 255/255]; %model 3

%Iterate through analysis for each trial we pick above
figure;
hold on;
for i=1:ntrials
    %Initiate Each Trial
    clear Data displacement load time
    filename=names(jj(i):kk(i));
    fileID=importdata(filename);
    j=jj(i)+numchar;
    k=kk(i)+numchar;
    jj(i+1)=j;
    kk(i+1)=k;
    Data=fileID.data;
    displacement=Data(:,1)-dat_offset(i);  %mm
    load=Data(:,2);  %N
    time=Data(:,3);  %s
    angles(i)=str2num(filename(16:19));
    Lo(i)=max(displacement);
    % Stress & Strain are not used but can be calculated as follows:
    %strain=(displacement./Lo(i)).*100; % %
    %stress=((load)./Area); %MPa
    plot((displacement-dat_offset(i))./plotnorm,load,'LineWidth',2','Color',colors(i,1,:));
    
    %Determine Peeling Model
    phi=angles(i);
    xmax=@(phi) 2.*w.*tand(phi./2);
    fun11=@(x) 2.*gamma.*(x+dat_offset(4)).*((cscd(phi./2)).^2).*(cotd(phi./2));
    fun12=@(x) 2.*gamma.*(xmax(phi)-(x+dat_offset(4))).*((cscd(phi./2)).^2).*(cotd(phi./2));
    xxl1=linspace(0,xmax(phi)./2,1000);
    xxl2=linspace(xmax(phi)./2,xmax(phi),1000);
    y1=(fun11(xxl1));
    y2=(fun12(xxl2));
    
    %Plot Peeling Model
    plot(xxl1./plotnorm,y1,'Color',colors(i,2,:),'LineWidth',2,'LineStyle','--');
    plot(xxl2./plotnorm,y2,'Color',colors(i,2,:),'LineWidth',2,'LineStyle','--');
    
    %Determine Peaks of tensile data
    text((xmax(phi)./2+text_offset(i))./plotnorm,((fun11(xmax(phi)./2)))+.2,[num2str(phi) '\circ'],'FontSize',14,'Color',colors(i,1,:));
    adjusttrials=[13 15 17 18 19 20 22 24]; 
    if any(ismember(adjusttrials,str2num(filename(12:14))))
        [pks, pkloc]=findpeaks(load,displacement);
        [pks2, pkloc2]=findpeaks(pks,pkloc);
        peaks=[pks2'; pkloc2']';
        sortpks=sortrows(peaks,1);
        peakloads(1:length(pks2),i)=sortpks(:,1);
        peakload(i)=peakloads(length(sortpks)-1,i);
        pklocs(i)=sortpks(length(sortpks)-1,2);
        adjustloads=[adjustloads max(pks)];
        adjustloc=[adjustloc sortpks(length(sortpks),2)-dat_offset(i)];
    else
        [pks, pkloc]=findpeaks(load,displacement);
        peaks=[pks'; pkloc']';
        sortpks=sortrows(peaks,1);
        peakload(i)=max(pks);
        pklocs(i)=sortpks(length(sortpks),2)-dat_offset(i);
    end
    xmaxs(i)=xmax(phi)./2;
    peelmax(i)=max(y1);
    
    %Determine Elastic Model
    xx=linspace(0,xmax(phi),1000);
    EFfun=@(x) -(1/2).*E.*t.*x.*tand(phi).*log(abs((2.*x.*cotd(phi./2)-beta.*tand(phi))./(-beta.*tand(phi))));
    EF=EFfun(xx);
    for n=1:(length(EF)-1)
        m=(EF(n+1)-EF(n))/(xx(n+1)-xx(n));
        b=EF(n)-m*xx(n);
        %linearize elastic model
        if m>slope
            xremain=linspace(xx(n),max(xx),(length(xx)-n+1));
            EF(n:length(xx))=xremain.*slope+b;
            break
        end
        linearpts(i)=xx(n);
    end
    
    %Plot Elastic Model
    plot(xx./plotnorm,EF,'Color',colors(i,2,:),'LineWidth',2,'LineStyle',':');
    
    %Determine Peeling Energy
    xmax=xmax(phi);
    Upeelinc=@(x) 2.*gamma.*((cscd(phi./2)).^2).*(1./2).*(x.^2).*(cotd(phi./2));
    Upeeldec=@(x) 2.*gamma.*((cscd(phi./2)).^2).*((1)./2).*((-1.*(xmax-x).^2)+2.*(xmax./2).^2).*(cotd(phi./2));
    xx=linspace(0,xmax,1000);
    xxl1=linspace(0,xmax./2,length(xx)./2);
    xxl2=linspace(xmax./2,xmax,length(xx)./2);
    xx=[xxl1 xxl2];
    Upeel=[Upeelinc(xxl1) Upeeldec(xxl2)];
    %Determine Elastic Energy
    Uelas=cumtrapz(xx,EF);
    
    %Determine energy inetersect
    subU=Uelas-Upeel;
    for a=1:length(subU)-1
        test=subU(a);
        if test>0
            break
        end
    end
    intersectUx(i)=xx(a);
    intersectUFelas(i)=EF(a);
end

%Plot Peeling Maxima
plot(xmaxs./plotnorm,peelmax,'Marker','.','Color',[180/255 0 0],'LineStyle','none','MarkerSize',30);

%Plot Max Load
scatter(pklocs./plotnorm,peakload,70,'filled','MarkerFaceAlpha',1/3,'MarkerFaceColor',[0 0 0]...
    ,'MarkerEdgeColor',[0 0 0],'MarkerEdgeAlpha',3/3);

%Plot Elastic Max 
scatter(adjustloc./plotnorm,adjustloads,70,'filled','MarkerFaceAlpha',3/3,'MarkerFaceColor',[255/255 139/255 130/255]...
    ,'MarkerEdgeColor',[180/255 68/255 54/255],'MarkerEdgeAlpha',3/3);

%Plot where Elastic Force meets alpha (orange square)
plot(linearpts./plotnorm,EFfun(linearpts),'Marker','s','MarkerFaceColor',[255/255 126/255 24/255],'LineStyle','none','MarkerSize',10);

%Plot Elastic Energy=Peeling Energy (yellow triange)
plot(intersectUx./plotnorm,intersectUFelas,'Marker','^','MarkerFaceColor',[255/255 193/255 37/255],'LineStyle','none','MarkerSize',10);
    
%Determine continuous Peeling Model
xmaxfun=@(phi) w.*tand(phi./2);
anglescontinuous=linspace(10,110,1000);
peelmax=zeros(1,length(anglescontinuous));
for j=1:length(anglescontinuous)
    phi=anglescontinuous(j);
    xmax=xmaxfun(phi);
    fun11=@(x) 2.*gamma.*(x).*((cscd(phi./2)).^2).*(cotd(phi./2));
    peelmax(j)=fun11(xmax);
end

%Plot  Continuous Peeling Model
plot(xmaxfun(anglescontinuous)./plotnorm,peelmax,'Color',[180/255 0 0],'LineWidth',2,'Marker','none');

%Formatting the plot
ylim([0 10])
xlim([0 1])
xlabel('x/w','FontSize',14);
ylabel('Load (N)','FontSize',14);
hold off;
ax=gca;
ax.FontSize=14;